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NOMENCLATURE 


X 

horizontal  (chordwise)  coordinate  of  physical  domain 

y 

vertical  coordinate  of  physical  domain 

r 

radial  coordinate  of  computational  domain 

e 

angular  coordinate  of  computational  domain 

J 

Jacobian  of  the  mapping 

K 

curvature  of  the  airfoil 

U 

x-component  of  velocity 

V 

y- component  of  velocity 

q 

the  velocity  vector  {u,v)  whose  magnitude  is  q 

<^ 

the  vorticity  vector,  V  x  q 

<!> 

velocity  potential 

c 

speed  of  sound 

p 

density 

p 

pressure 

a 

angle  of  attack 

r 

circulation 

Cd 

coefficient  of  drag 

CwD 

coefficient  of  wave  drag 

Cl 

coefficient  of  lift 

Cp 

coefficient  of  pressure 

Cv 

specific  heat,  constant  volume 
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internal  energy 

E 

total  energy,  q^ /2-\-  e 
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specific  enthalpy,  e  +  p/ p 
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total  enthalpy,  5^/2  +  i 

S 

entropy 

T 

temperature 

7 

adiabatic  exponent 

M 

Mach  number,  q/c 

e* 

momentum  thickness 

8 

displacement  thickness 

H* 

shape  factor 
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CHAPTER  I 


INTRODUCTION 


It  is  well  known  that  the  most  efficient  cruising  speed  of  modern  commercial 
aircraft  lies  in  a  range  of  free  stream  Mach  numbers  Moo  near  .7,  where  the  fluid 
flow  past  the  body  may  be  characterized  as  transonic  or  supercritical.  The  flight 
efficiency  is  measured  by  the  dimensionless  parameter  Moo  Cl/Cd,  where  Cl  and 
Cd  are  the  coefficients  of  lift  and  drag,  respectively.  In  this  range  of  Mach  numbers 
a  region  of  supersonic  flow  forms  above  the  wing,  which  is  terminated  by  a  shock 
wave  across  which  the  pressure,  density,  and  velocity  of  the  fluid  change  rapidly. 
A  rapid  rise  in  the  drag  due  to  these  shock  waves,  which  reduces  the  overcdl  flight 
efficiency,  occurs  in  the  transonic  range.  The  goal  of  the  aeronautical  engineer  is 
to  design  supercritical  airfoils  that  will  postpone  the  onset  of  strong  shocks  until  a 
higher  cruising  speed  can  be  attained. 

The  transonic  problem  originated  in  the  late  1940's  when  wind  tunnel  exper- 
iments seemed  to  indicate  that  shockless  supercritical  flows  do  not  exist.  Math- 
ematically this  suggested  that  the  boundary  value  problem  for  transonic  flow  is 


ill-posed  in  the  sense  that  no  smooth  solution  to  the  governing  equations  with  ar- 
bitrary boundary  data  exists.  In  1956,  Morawetz  [24]  proved  that  the  problem  is 
indeed  ill-posed.  However,  it  was  not  until  about  ten  years  later  that  the  mathe- 
matics community  began  to  accept  the  idea  that  there  should  always  exist  a  unique 
weak  solution  to  the  transonic  boundary  value  problem  if  an  appropriate  entropy 
inequality  is  imposed  [10].  In  the  mean  time,  the  problem  was  left  largely  to  the 
experimentalists  [29,33],  whose  successes  eventually  led  the  mathematicians  to  re- 
consider the  role  of  weak  solutions. 

With  the  development  of  faster  and  more  efficient  computers,  it  became  possi- 
ble to  formulate  more  realistic  numerical  models  of  transonic  flow.  At  first  research 
was  limited  to  computing  solutions  of  the  transonic  small  disturbance  equation,  in 
which  it  is  assumed  that  the  flow  is  isentropic,  irrotational,  and  only  slightly  per- 
turbed from  its  free  stream.  By  the  late  1960's,  Bauer,  Garabedian  and  Korn  [3,5] 
developed  the  method  of  complex  characteristics  for  the  design  of  shockless  airfoils, 
in  which  the  full  potential  equation  is  treated  in  the  hodograph  plane,  rendering  it 
linear  in  the  transformed  variables.  Comparable  work  had  been  published  by  Nieuw- 
land  [27]  at  the  NLR  in  Amsterdam.  In  1971,  Murman  and  Cole  [25]  introduced 
a  method  of  shock  capturing  to  solve  the  small  disturbance  equation  by  employing 
upwind  differencing  in  the  supersonic  zone  while  retaining  central  differencing  in 
the  subsonic  zone. 

These  ideas  led  to  codes  which  could  analyze  airfoils  in  either  two  or  three 


dimensions  by  solving  the  full  potential  equation.  In  particular,  the  BGK  code 
written  by  Bauer,  Garabedian,  and  Korn  gained  wide  acceptance  in  the  aeronautics 
community  and  elided  in  the  design  of  commercial  supercritical  airfoils  [4].  Figure  1 
displays  the  results  of  a  typical  analysis  run  using  the  BGK  code.  Presented  are  both 
computed  and  experimental  pressure  profiles  and  important  flow  quantities,  along 
with  the  sonic  line  bounding  the  supersonic  zone  and  a  long  hash-mark  on  the  Cp 
axis  denoting  the  critical  pressure.  This  particular  example  may  be  characterized 
as  nearly  shockless;  that  is,  the  airfoil  geometry  was  derived  through  the  method 
of  complex  characteristics  described  above,  and  the  flow  parameters  M  and  Cl  are 
close  to  their  design  values.  Reference  [4]  presents  the  results  of  numerical  tests 
performed  on  a  wide  variety  of  airfoils,  in  which  comparisons  with  experimental 
data  verify  the  accuracy  of  the  method. 

As  good  as  the  potenticd  model  is,  some  still  question  the  isentropic  assumption. 
It  has  been  argued  that  the  jump  in  entropy  across  a  shock  may  not  be  neglected 
if  its  strength  and  location  are  to  be  computed  accurately.  The  Euler  equations 
constitute  a  more  mathematically  correct  model,  and  when  written  in  conservative 
form,  they  correctly  treat  the  shock  jumps.  By  the  early  1980's,  several  codes  were 
introduced  to  compute  the  solution  to  the  Euler  equations.  In  particular,  the  two- 
dimensional  code  FL052  written  by  Jameson  [14]  was  highly  successful.  The  scheme 
is  relatively  easy  to  implement  because  the  equations  are  written  in  conservative 
form,  all  spatial  derivatives  are  centrally  differenced,  and  the  solution  is  iterated  in 


time  through  a  standard  Runge-Kutta  fourth  order  scheme.  Its  artificial  viscosity 
terms  are  a  blend  of  second  and  fourth  order  derivatives,  closely  resembling  the 
Laplacian  and  the  biharmonic  operator,  and  the  robust  mechanism  by  which  they 
are  switched  on  and  off  is  the  key  to  the  success  of  the  method. 

To  evaluate  the  merits  of  the  Euler  formulation,  we  have  incorporated  the 
Jameson  scheme  into  the  BGK  potential  code  to  study  the  two  models  in  a  common 
framework.  It  is  thus  possible  to  compare  the  solutions  of  both  the  nonconservative 
and  conservative  forms  of  the  potential  equation  directly  with  the  solution  of  the 
Euler  equations.  We  have  studied  how  the  three  formulations  treat  the  problem  of 
flow  past  an  airfoil  at  realistic  operating  conditions  where  the  wave  drag  is  not  large 
enough  to  cause  the  flow  to  separate.  In  this  range  our  comparison  has  shown  that 
the  potential  and  the  Euler  models  are  more  or  less  equivalent. 

In  fact,  for  flow  simulations  with  a  boundary  layer  correction,  we  have  found 
that  the  solution  of  the  nonconservative  potential  equation  actually  matches  exper- 
imental data  better  than  the  solution  of  the  Euler  equations.  The  reason  for  this 
surprising  result  stems  from  the  treatment  of  shock  wave-boundary  layer  interac- 
tion by  the  nonconservative  method.  This  had  been  noted  for  some  time  in  regard 
to  comparisons  between  the  nonconservative  and  conservative  potential  equation, 
where  it  had  been  observed  that  the  former  method  provides  a  better  estimate  of 
shock  strength  and  location  [4,22].  It  might  be  argued  that  the  isentropic  assump- 
tion is  the  cause  of  the  inaccuracy  of  the  conservative  potential  solution,  and  that 
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it  is  necessary  to  solve  the  Euler  equations  in  conservative  form  to  arrive  at  the 
proper  solution.  However,  we  will  show  that  the  cause  of  the  discrepancy  is  not  the 
isentropic  assumption,  but  rather  shock  wave-boundary  layer  interaction. 

Obviously,  it  is  crucial  to  this  study  that  we  be  able  to  calculate  the  wave 
drag  for  the  Euler  solution  accurately.  For  computations  with  a  boundary  layer 
correction,  the  standard  pressure  integration  around  the  airfoil  is  unreliable  due  to 
the  difficulty  in  computing  the  flow  variables  in  the  wake.  As  an  alternative,  we  will 
show  that  the  Euler  equations  modified  by  Jameson's  artificial  viscosity  terms  can 
be  combined  to  form  an  equation  for  the  entropy  whose  viscous  terms  provide  us 
with  a  means  of  calculating  the  wave  drag.  Instead  of  differencing  uncertain  values 
of  pressure,  we  integrate  a  positive  definite  quantity  which  depends  on  the  viscous 
coefficients  and  the  gradients  of  the  flow  variables  over  a  region  surrounding  the 
shock,  and  thereby  we  avoid  the  trouble  at  the  nose  and  in  the  wake. 

We  include  in  the  final  chapter  a  discussion  of  how  the  Euler  scheme  handles 
vorticity.  This  is  not  only  interesting  in  its  own  right,  but  it  also  provides  us 
with  insight  into  the  biharmonic  terms  of  the  artificial  viscosity.  Even  though 
the  Euler  equations  admit  an  infinite  number  of  rotational  solutions  in  which  closed 
streamlines  may  occur,  the  dissipative  nature  of  the  biharmonic  terms  in  most  cases 
forces  the  scheme  to  converge  to  the  irrotational  or  nearly  irrotational  solution.  This 
is  certainly  the  case  for  airfoils  and  other  streamlined  bodies.  Also,  we  will  show  that 
the  Euler  scheme  does  not  automatically  satisfy  the  Kutta- Joukowski  condition,  but 
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rather  relies  on  the  biharmonic  terms  to  diffuse  vortices  which  are  shed  from  the 
sharp  traihng  edge  of  an  airfoil.  This  realization  explains  some  of  the  difficulties 
encountered  when  computing  the  flow  variables  in  the  wake. 

Beyond  apphcations  to  cdrfoils,  we  have  found  the  Euler  scheme  to  be  an  ex- 
tremely useful  tool  in  the  study  of  vortex  dynamics.  With  no  changes  to  the  fun- 
damental algorithm,  we  are  able  to  compute  rotational  flows  over  a  wide  variety 
of  blunt  bodies.  By  initializing  the  scheme  with  known  incompressible  rotational 
flows,  we  are  able  to  track  vortices  rather  well  —  both  upstream  and  downstream, 
depending  on  their  strength.  We  are  able  to  compute  the  flow  around  a  plate  which 
is  perpendicular  to  the  free  stream.  By  increasing  the  free  stream  Mach  number, 
we  are  able  to  study  the  interaction  of  shocks  with  regions  of  high  vorticity.  We  are 
even  able  to  simulate  such  complex  phenomena  as  the  von  Karman  vortex  street 
and  the  Helmholtz  instability  on  a  non-uniform  and  relatively  crude  mesh. 


CHAPTER  II 
MODELS  OF  TRANSONIC  FLOW 

1.  Equations  of  Gas  Dynamics 

The  two-dimensional  motion  of  an  ideal,  inviscid  fluid  is  governed  by  the  Euler 
equations 

(2.1a)  pt  +  {pu)x  +  ipv)y  =  0, 

(2.1b)  {pu)t  +  {pu^  +  p)^  +  {puv)y  =  0, 

(2.1c)  {pv)t  +  (puv)^  +  (pv^  +  p)y  ^  0, 

(2.1d)  {pE)t  +  (puH)^  +  {pvH)y  =  0, 

where  E  =  q"^ /2  +  e,  and  H  =  E  +  p/p.  The  fluid  we  are  considering  is  air  at 
moderate  temperatures,  whose  thermodynamic  state  is  completely  defined  by  the 
quantities  p,  p,  e,  5,  and  T.  We  choose  to  express  the  last  three  quantities  in  terms 
of  the  first  two  by  taking  the  gas  to  be  polytropic;  that  is,  its  interned  energy  e  is 
assumed  to  be  directly  proportional  to  its  temperature  T.  This  leads  to  an  equation 
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of  state 


(2.2)  p^A{S)p\ 


where  7  =  1.4  for  air  and  the  function  A{S)  is  given  as 


(2.3)  ^(5)  =  ^exp^*^      ^"^ 


Poo  \  Cyj 

The  quantities  e  and  T  may  be  expressed  in  terms  of  p  and  p  by 

(2.4)  e  =  c„r=-^^. 

7-  1  p 

Equations  (2.1)  along  with  the  polytropic  assumption  constitute  what  will  be 
referred  to  as  the  Euler  model  of  fluid  flow.  The  potential  model,  in  addition, 
takes  the  flow  to  be  in  steady  state,  removing  derivatives  with  respect  to  time 
t  in  equations  (2.1);  and  it  assumes  that  the  entropy  S  is  constant  throughout 
the  flow  domain,  eliminating  equation  (2. Id)  and  rendering  the  function  A{S)  a 
constant.  Because  the  flow  is  uniform  at  infinity,  Kelvin's  theorem  asserts  that  the 
flow  is  irrotational,  and  thus  there  exists  a  velocity  potential  (f)  whose  gradient  is  the 
velocity  vector  q.  From  these  assumptions  we  obtain  a  single  second  order  equation 
for  the  velocity  potential  cf) 

(2.5)  (C^    -   <i>l)  <^„    -   2  <i>.cl>y<t>,y   +  (C2    -   Ci>l)4>yy    =   0. 

This  equation,  which  will  be  referred  to  as  the  nonconservative  potential  equation, 
is  elliptic  when  M^  <  1  and  hyperboHc  when  KP  >  1.  Equation  (2.5)  may  also  be 
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written  in  conservative  form  as 

(2.6)  (p</>x)x  +  (p<A«)y  =  0. 

This  will  be  referred  to  as  the  fully  conservative  potential  equation.  Both  forms  of 
the  potential  equation  make  use  of  Bernoulli's  law 

(2.7)  \q^  +  -:L-l^H, 

2  7-  Ip 

where  Hy  the  total  enthalpy,  is  constant  throughout  the  flow.  Bernoulli's  law  is  a 
consequence  of  equations  (2.1)  in  steady  state,  the  uniformity  of  the  flow  at  infinity, 
and  the  conservation  of  H  across  a  shock.  This  last  fact  will  be  discussed  in  Section 


2.  Boundary  Conditions 

Both  models,  the  Euler  and  the  potential,  employ  a  no-flux  boundary  condition 
at  the  airfoil  surface  for  an  inviscid  calculation,  or  at  the  displaced  airfoil  surface 
for  a  viscous  calculation  employing  a  boundary  layer  correction.  For  the  potential 
equation,  this  is  equivalent  to  the  Neumann  boundary  condition 

(2,8)  1^=0. 

For  the  Euler  equations,  the  no-flux  condition  asserts  that  the  normal  component 
of  the  velocity  is  zero  at  the  airfoil  surface.  Consequences  of  this  simple  condition 
for  the  Euler  model  will  be  described  in  Chapter  III,  Section  4. 
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In  addition  to  the  no-flux  boundary  condition,  the  solution  must  satisfy  the 
Kutta-Joukowski  condition,  which  asserts  that  the  flow  velocity  remains  finite  at 
the  traihng  edge  of  the  airfoil.  For  the  potential  model  the  circulation  around  the 
airfoil  F  is  given  by  the  jump  in  (j)  across  the  traihng  streamline  and  is  easily  adjusted 
to  move  the  separation  point  to  the  tail.  The  Euler  model  has  no  explicit  mechanism 
to  enforce  the  Kutta  condition  [31],  relying  instead  on  its  artificial  viscosity  terms 
to  suppress  the  flow  from  rounding  the  traihng  edge.  This  issue  will  be  taken  up 
again  in  Chapter  VI. 

The  far  field  boundary  condition  for  the  potential  model  is  also  easily  imple- 
mented. It  can  be  shown  [19]  that  </>  has  an  asymptotic  expansion 


r 

(2.9)  (j)  ~  qooR  cos(0  -  q)  +  —  tan-^  (l3  tan(0  -  a)) 

ZTT 


as  R'^  =  x^  +  y"^  ^  oo;  where  ^"^  —  1  -  M^,  Q  =  ta.n~^{y/x)  and  a  is  the  angle 
of  attack.  In  order  to  solve  the  potential  equations  numerically,  the  physical  flow 
domain  is  mapped  conformally  into  the  unit  circle,  with  the  point  R  =  oo  sent 
to  the  origin.  The  boundary  condition  (2.9)  may  then  be  easily  imposed  once  4> 
has  been  modified  to  remove  the  singularity  associated  with  the  first  term  in  the 
expansion  [11].  Following  each  potential  iteration  the  circulation  F  is  updated  to 
coincide  with  the  value  obtained  through  the  enforcement  of  the  Kutta-Joukowski 
condition  on  the  airfoil. 

The  flow  domain  for  the  Euler  model,  however,  does  not  extend  to  infinity,  and 
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thus  at  the  far  field  boundary  we  must  implement  a  radiation  condition  based  on 
the  number  of  incoming  and  outgoing  characteristics.  In  Chapter  III,  we  discuss 
this  matter  in  greater  detail. 

3.   Shock  Conditions 

For  subsonic  flow  the  potential  and  the  Euler  formulations  of  the  equations  of 
motion  are  equivalent,  but  for  transonic  flow  the  two  models  yield  diff'erent  shock 
structures.  Denoting  by  N  and  L  the  normal  and  tangential  components  of  the  ve- 
locity vector  q  at  a  stationary  shock  front,  we  have  the  Rankine-Hugoniot  relations 
corresponding  to  the  full  set  of  Euler  equations: 

(2.10a)  [pN]  =  0, 

(2.106)  [pN''+p]^0, 

(2.10c)  [L]  =  0, 

(2.10cf)  [H]  =  0. 

Here,  the  bracket  [  •  ]  represents  a  jump  in  the  enclosed  quantity  across  a  shock.  To 
these  relations  we  add  a  fifth  condition  which  is  a  consequence  of  the  viscous  effects 
which  may  not  be  neglected  in  the  presence  of  steep  gradients  in  the  flow  variables. 
It  is  referred  to  as  the  entropy  inequality 

(2.11)  [S]  >  0, 

and  it  asserts  that  the  entropy  of  the  fluid  increases  as  it  passes  through  a  shock 
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wave.  It  can  be  shown  [7]  that  the  jump  in  entropy  is  of  third  order  in  the  strength 
of  the  shock,  where  the  strength  is  measured  by  the  jump  in  the  pressure. 

The  isentropic  assumption  of  the  potential  model  does  not  allow  for  a  jump 
in  entropy  across  a  shock,  and  this  forces  us  to  give  up  one  of  the  first  four  shock 
conditions  (2.10),  specifically  the  conservation  of  normal  momentum  (2.10b).  The 
role  of  the  entropy  in  the  Euler  model  is  replaced  by  the  component  of  normal 
momentum  in  the  potential  model;  in  other  words,  the  condition 

(2.12)  [/>iV2+p]>0 

forces  the  shocks  to  be  compressive  [16].  Thus,  in  the  potential  model,  shocks  create 
a  momentum  excess  which,  in  the  case  of  flow  past  a  wing,  serves  as  a  retarding 
force  on  the  body.  Garabedian  [9]  makes  use  of  this  phenomenon  in  calculating  the 
wave  drag.  This  subject  will  be  taken  up  in  Chapter  IV. 
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CHAPTER  III 


THE  BGKM  COMPUTER  CODE 


1.  Introduction 

In  this  chapter  we  present  some  of  the  technical  details  associated  with  the 
BGKM  code.  We  are  primarily  concerned  with  the  implementation  of  an  Euler 
solver  as  an  option  in  the  original  potential  code.  This  scheme  is  based  on  the 
work  of  Jameson,  Schmidt,  and  Turkel  [13].  Details  of  the  potential  routines  may 
be  found  in  Reference  [4].  Further  documentation  concerning  the  operation  of  the 
code  may  be  found  in  the  listing  itself,  which  is  available  on  request. 

We  should  mention  that  the  Euler  scheme  which  we  have  implemented  does 
not  include  some  of  the  features  which  have  been  suggested  by  the  original  authors 
to  speed  up  the  convergence  time.  These  include  enthalpy  damping  and  residual 
averaging  [14].  For  our  study  we  do  not  necessarily  require  the  fastest  Euler  solver, 
and  the  additional  machinery  required  to  speed  up  the  code  might  obscure  the  role 
played  by  the  artificial  viscosity  terms  which  are  added  to  the  Euler  equations. 
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w  =  .    .  ,  _  ,  pu^+P 


These  terms  are  of  vital  interest,  as  we  shall  see  in  Chapters  IV  through  VI. 

2.   Modified  Euler  Equations 

The  Euler  equations  in  conservative  form  may  be  written  as 

(3.1)  wt  +  f^+gj,  =  0, 

where  x  and  y  are  Cartesian  coordinates  and 

pu      \  /      pv      \ 

I      puv 

,  g    =  9 

puv  I   pv    +  p 

puH    J  \    pvH    I 

In  order  to  solve  equations  (3.1)  numerically,  we  first  map  the  region  exterior  to  the 

airfoil  into  the  unit  circle,  where  the  equations  may  be  written  again  in  conservative 

form  in  terms  of  polar  coordinates  r  and  6.  Then  artificial  viscosity  terms  are  added 

to  the  right  hand  side,  resulting  in  a  new  set  of  equations  to  be  referred  to  as  the 

modified  Euler  equations 

(3.2)  Wt  +  Fe  +  G,  =  V, 

where  W  =  J  w,  F  =  t/r  f  —  x^  g,  G  =  a;^  g  —  y©  f ,  J  =  xgyr  —  Xrl/o,  and  the 
artificial  viscosity  vector  V  will  be  described  in  detail  presently.  Derivatives  with 
respect  to  r  and  6  are  approximated  using  central  differences,  and  the  flow  vector  W 
is  advanced  in  time  through  a  standard  Runge-Kutta  fourth  order  scheme.  Details 
will  be  presented  in  Section  3. 
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Unlike  the  finite-volume  formulation  of  the  Euler  equations  solved  by  Jameson, 
Schmidt  and  Turkel  [13],  the  equations  (3.2)  are  written  in  finite  difference  form 
to  conform  with  the  established  conventions  of  the  original  BGK  potential  code  of 
Bauer,  Garabedian  and  Korn.  MacCormack  [20]  demonstrates  the  equivalence  of 
the  two  formulations  at  interior  mesh  points,  noting  that  the  equivalence  is  not  valid 
at  boundaries.  Finite  difference  methods  place  mesh  points  on  boundaries,  whereas 
finite  volume  methods  place  surface  elements  along  boundaries.  It  should  be  noted 
that  the  mesh  point  which  falls  on  the  sharp  trailing  edge  must  be  handled  with 
special  care.  Treatment  of  the  boundaries  will  be  discussed  in  Section  4. 

The  artificial  viscosity  vector  which  is  added  to  the  right  hand  side  of  equation 
(3.2)  is  a  combination  of  second  and  fourth  order  finite  difference  terms  in  both  the 
angular  and  radial  directions,  which  are  indexed  by  the  integers  /  and  j ,  respectively. 
The  radial  terms  take  the  form 

(3.3)        v^")  =  d;  ( u^l,  Dt )  w'  -  d;  ( u\';l,  Dt  d;  Dt )  w', 

where  D^  and  D~  denote  forward  and  backward  finite  differences,  respectively.  The 
angular  terms  are  similarly  defined.  Here  w'  is  the  same  as  w  except  that  pE  is 
replaced  by  pH  in  the  fourth  row  to  ensure  that  the  total  enthalpy  H  is  conserved 
even  after  the  addition  of  the  artificial  viscosity  terms.  It  should  be  noted  that 
the  artificial  viscosity  operator  may  act  on  either  pE  or  pH  in  the  fourth  Euler 
equation  —  we  shall  see  that  the  entropy  inequality  is  satisfied  in  both  cases.  We 
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have  chosen  to  act  on  pH  so  that  output  values  of  H  provide  a  convenient  check  on 
the  convergence  and  accuracy  of  the  method. 

The  coefficient  v^"^^  of  the  second  order  difference  terms  of  (3.3)  is  defined  as 

(3-4)  t.[^/^,  =(J/AO,j+i   e,,,+  i, 

where 


(J/A<),,,+  i   =  \{{Jl^t)i,i  +  (J/AO/.i+i), 


t,  •  I  1    —   K2      max 


P/,i-i-i  -2pi,i  +  pi,i-i 


P/,i+i  +  2pi,i  +  p/,i-i 

A<  is  the  locaJ  time  step  determined  by  calculating  the  CFL  condition  at  each 
point  in  the  flow  domain.  The  key  component  of  the  the  viscosity  coefficient  u^"^' 
is  the  term  involving  the  second  derivative  of  pressure  which  essentially  signals  the 
formation  of  a  shock.  Consequently,  the  second  order  difference  terms  become  of 
first  order  in  the  grid  spacing  near  a  shock  and  of  third  order  away  from  a  shock. 
As  will  be  shown  in  the  next  chapter,  this  ensures  that  the  method  satisfies  the 
entropy  inequality. 

The  coefficient  u^*^  of  the  fourth  order  difference  terms  of  (3.3)  is  given  as 

(3.5)  i/fj^i    =  (-//AOij+i   max(0,K4-e,j+i). 

The  fourth  order  terms  are  important  in  regions  where  the  flow  is  relatively  smooth, 
where  they  are  of  third  order  in  the  grid  spacing;  and  are  switched  off  near  a  shock  to 
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prevent  overshoots.  In  Chapter  VI  we  will  show  that  these  terms  not  only  stabiUze 
the  scheme  in  the  far  field,  but  that  they  also  are  needed  in  the  wake  to  enforce  the 
Kutta-Joukowski  condition  at  the  traihng  edge. 

The  parameters  K2  and  K4  of  (3.4)  and  (3.5)  are  adapted  to  the  flow.  Typical 
values  for  an  airfoil  are  .5  and  .015,  respectively.  Higher  values  are  used  for  a  blunt 
body  such  as  a  cylinder  or  stalled  wing.  In  Chapter  VI  we  will  investigate  the  effects 
on  flow  computations  caused  by  varying  these  two  parameters. 

3.  Euler  Time-Stepping  Scheme 

The  local  time  step  At  is  determined  each  iteration  and  at  each  point  in  the 
domain  according  to  the  CFL  condition 

(3.6)  A^      <      "^ ^, 

Ax  ^  Ay  ^  ^   \  Ax  ^  Ay  J 

This  CFL  bound,  suggested  by  Jameson  [14],  is  based  on  a  one-dimensional  model 
equation  similar  to  the  linearized  Burgers  equation.  Because  we  are  usually  inter- 
ested in  the  steady  state  solution  of  the  modified  Euler  equations,  the  solution  is 
advanced  in  non-physical  time  to  speed  up  the  method.  This  idea  will  be  of  concern 
to  us  in  Chapter  VI  where  time-dependent  solutions  are  considered. 

All  spaticd  derivatives  of  equation  (3.2)  are  centrally  differenced,  and  the  so- 
lution is  iterated  in  time  through  a  standard  fourth  order  Runge-Kutta  stepping 
scheme.  Note  in  the  following  that  the  artificial  viscosity  vector  V(w(''^)  is  com- 
puted only  once  during  each  iteration  to  speed  up  the  procedure,  w^"^  is  updated 
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according  to  the  scheme 


^(0)   ^  ^(n) 


w(i)=w(°)-^P(w('')) 
(3.7)  w(2)=w(«)-— P(w(i)) 


Ai 


^(n+i)  ^  ^(0)  _  ^  (|p(w(o))  +  2P(w(^))  +  2P(w(2>)  +  P(w(3))), 


w 


here 


(3.8)  P(w(''>)  =  -  ('f(w('=))<,  +  G(w(*))^  -  V(w(°)))  . 

References  [6,8]  consider  other  forms  of  the  Runge-Kutta  time-stepping  scheme, 
especially  in  conjunction  with  adaptive  and  multigrid  methods.  The  Euler  solver 
added  to  the  BGKM  code  does  not  yet  have  a  multigrid  capabihty,  but  normally  the 
potential  solvers  provide  a  good  enough  initialization  to  enable  the  Euler  solution 
to  be  computed  rapidly  on  a  fine  mesh.  The  potential  solvers,  however,  do  employ 
a  simplified  multigrid  scheme  by  solving  first  on  a  crude  grid,  and  then  on  a  grid 
twice  as  fine.  In  Section  4  we  will  describe  the  procedure  by  which  the  Euler  solver 
may  be  initialized  with  the  potential  solution. 
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4.  Euler  Boundary  Conditions 

At  the  airfoil  surface,  we  must  enforce  the  no-flux  condition, 

(3.9)  qn  =  0. 

Since  the  transformed  Euler  equations  (3.2)  are  written  in  conservative  form,  only 
the  pressure  is  needed  to  compute  the  vector  G  at  the  solid  boundary  because  all 
terms  involving  the  normal  flux  drop  out.  The  fact  that  only  the  pressure  need  be 
specified  agrees  with  the  fact  that  only  one  characteristic  carries  information  to  the 
impermeable  boundary  from  the  interior  of  the  flow.  The  pressure  on  the  boundary 
may  be  linearly  extrapolated  from  the  value  of  the  pressure  at  the  first  interior  grid 
point,  using  the  normal  derivative  of  the  pressure  at  the  boundary 

(3.10)  ^  =  ^pq\ 

where  k  is  the  curvature  of  the  airfoil.  Equation  (3.10)  is  derived  from  the  normal 
momentum  equation  at  the  boundary 

(3.11)  n-(pqt  +  pq-Vq+ Vp)  =  0, 
which  may  be  rewritten  due  to  (3.9)  as 

(3.12)  Vp.n  =  /)q.(q.V)n. 

This  extrapolation  of  pressure  is  not  apphcable  at  the  sharp  trailing  edge  due  to  the 
infinite  curvature  at  this  point.  However,  since  the  Jacobian  of  the  mapping  is  zero 
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there,  no  contribution  is  made  to  the  momentum  flux  in  the  direction  of  the  traihng 
streamline,  and  thus  the  pressure  need  not  be  specified  at  this  point.  (The  vector 
G  of  equation  (3.2)  is  identically  zero  there.)  Values  of  q  and  p  at  the  boundary, 
needed  only  for  a  boundary  layer  correction,  are  computed  by  assuming  that  the 
flow  is  radially  isentropic,  i.e.  that  there  is  no  diff^erence  in  the  entropy  between 
a  boundary  point  and  the  point  just  inside  the  flow  domain.  This  assumption  is 
reasonable  because  large  gradients  in  the  entropy  occur  only  at  shocks  and  are 
directed  parallel  to  the  boundary. 

At  each  interior  grid  point,  the  spatial  derivatives  of  equation  (3.2)  are  com- 
puted from  the  values  of  the  flow  vector  w  one  grid  point  away  in  each  coordinate 
direction.  The  artificial  viscosity  terms  (3.3),  however,  require  a  knowledge  of  w 
two  grid  points  away  in  each  direction.  Thus,  at  grid  points  adjacent  to  the  bound- 
ary points,  attention  must  be  paid  to  the  evaluation  of  the  viscous  terms.  Near  the 
airfoil  boundary,  radial  diff^erences  in  the  vector  (3.3)  are  simply  set  to  zero.  This 
change  does  not  affect  the  computation  along  the  top  or  bottom  of  the  airfoil  where 
the  steep  gradients  are  in  the  angular  direction,  but  it  can  be  quite  destabilizing 
at  the  tail.  The  stalled  wing  simulations  of  Chapter  VI  point  out  the  reasons  for 
such  difficulties.  Unlike  potential  flow,  an  Euler  flow  with  a  shock  will  exhibit  a 
discontinuity  in  velocity  at  the  trailing  edge  forming  a  slip  line.  In  a  real  flow,  this 
slip  line  is  unstable  and  roUs  up  to  form  the  turbulent  wake.  This  is  difficult  to  treat 
numerically,  and  is  made  more  difficult  due  to  the  trouble  computing  the  artificial 
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viscosity  vector  in  the  wake.  In  Chapter  VI  we  will  take  up  this  issue  again  when 
we  consider  the  implicit  Kutta-Joukowski  condition. 

Far  field  boundary  conditions  are  more  difficult  to  implement.  At  the  outer 
boundary  the  Euler  equations,  under  a  rotation  of  coordinates  into  the  radial  and 
angular  directions,  yield  four  Riemann  invariants  in  the  normal  (radial)  direction 
if  derivatives  in  the  angular  direction  are  neglected.  Let  A^  denote  the  outward 
normal  component  of  the  velocity  vector  q  at  the  outer  boundary  and  L  the  lateral 

component.  Then  we  have  the  following  characteristic  equations: 

dS  dx 

^=0  aJong     -  =  Ar, 

dL  dx 

(3.X3)  ,,  ,t  ""     *^''' 

Here  x  denotes  the  outward  normal  component  of  the  local  orthogonaJ  coordinate 
system.  Where  the  fiow  is  inward,  that  is,  coming  from  infinity  into  the  flow  do- 
main, there  are  three  incoming  and  one  outgoing  characteristics,  and  where  the 
flow  is  outward,  there  are  three  outgoing  and  one  incoming  characteristics.  The 
Riemann  invariants  corresponding  to  incoming  characteristics  take  on  free  stream 
values,  whereas  those  corresponding  to  outgoing  characteristics  take  on  values  ex- 
trapolated from  the  flow  domain.  Each  extrapolated  quantity  is  computed  as  a 
weighted  average  of  its  value  at  the  outer  most  grid  point  and  at  the  next  point 
inside  the  flow  domain.  The  weighting  reflects  the  characteristic  speed  of  each  Rie- 
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mann  invariant.  We  are  free  to  assign  the  pressure  and  the  velocity  of  the  free  stream 
the  value  of  one.  Then  the  density  may  be  determined  from  the  thermodynamic 
relations 


(3.14)  cl    -   ^°°    -  ^^°° 


'oo 


Ml 


and  the  components  of  velocity  N  and  L  may  be  determined  from  the  angle  of 
attack. 

The  procedure  described  above  essentieilly  asserts  a  radiative  condition  at  the 
outer  boundary,  but  does  not  provide  the  equivalent  of  the  far  field  boundary  con- 
dition of  the  potential  model  (2.9).  There  is  no  assertion  of  a  circulation  at  infinity 
—  it  is  assumed  that  the  circulation  wrill  be  estabHshed  automatically  in  the  far 
field  if  the  outer  boundary  is  sufficiently  far  from  the  airfoil.  Also,  the  fact  that 
the  free  stream  values  are  used  by  the  fourth  order  difference  terms  of  the  artificial 
viscosity  vector  (3.3)  necessitates  that  the  outer  boundary  be  far  enough  away  from 
the  airfoil  so  that  the  flow  may  be  assumed  to  be  unperturbed  from  its  free  stream 
state  there.  In  the  far  field,  the  fourth  order  differences  are  crucial  for  maintaining 
numerical  stability,  and  differences  in  the  radial  direction  cannot  simply  be  switched 
off,  but  rather  must  be  formed  using  free  stream  values  of  w.  For  these  reasons, 
Thomas  and  Salas  [32]  suggest  that  the  outer  boundary  should  be  at  least  50  chord 
lengths  from  the  airfoil,  although  a  simple  numerical  test  suggests  that  25  chord 
lengths  is  sufficiently  far.  We  computed  hfting  solutions  of  the  Euler  equations  for 
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meshes  of  various  radicil  extents,  and  found  that  no  significant  loss  in  lift  occurred 
for  meshes  extending  at  least  25  chord  lengths  from  the  airfoil.  If  the  radial  extent 
is  lessened  to  10  or  15  chord  lengths,  significant  losses  in  the  computed  lift  are  ob- 
served. Reference  [32]  presents  a  procedure  in  which  the  leading  order  term  of  the 
small  disturbance  equation  is  incorporated  into  the  far  field  boundary  conditions 
for  meshes  extending  as  little  as  5  chord  lengths  from  the  body. 

5.  Code  Operation 

The  BGKM  computer  code  can  both  design  and  analyze  supercritical  airfoils. 
In  the  design  mode,  the  code  takes  as  input  a  desired  pressure  distribution  around 
an  airfoil  and  proceeds  to  solve  the  free-boundary  problem  which  produces  the 
coordinates  of  the  airfoil  associated  with  such  a  pressure  profile.  Reference  [21] 
contains  a  detailed  description  of  this  procedure.  The  problem  we  are  considering 
is  one  of  analysis,  rather  than  design.  In  fact,  the  design  mode  of  the  code  makes 
no  use  of  the  Euler  solving  apparatus  described  above. 

It  is  extremely  informative  to  compare  the  models  of  transonic  flow  by  plotting 
the  pressure  profile  of  each  on  a  single  graph.  The  BGKM  code  solves  the  noncon- 
servative  (NC)  potential,  the  fuUy  conservative  (FC)  potential,  and  the  modified 
Euler  equations  using  previously  computed  solutions  to  initialize  the  subsequent 
computations.  A  typical  analysis  run  involving  all  the  methods  proceeds  as  follows: 
1.  Given  airfoil  coordinates,  the  modulus  of  the  derivative  of  the  conformal  map- 
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ping  function  is  computed  for  use  by  the  potential  equations  written  in  the 
transformed  coordinates  r  and  9. 

2.  Due  to  its  use  of  a  fast  Poisson  solver,  the  solution  of  the  NC  potential  equation 
is  quickly  computed,  using  the  incompressible  solution  as  an  initial  guess. 

3.  The  FC  potential  solution  is  computed,  using  the  NC  potential  solution  to 
initialize.  Both  the  NC  and  FC  potential  solvers  employ  a  simple  multigrid 
technique,  solving  the  potential  equation  first  on  a  crude  grid,  then  on  one 
twice  as  fine. 

4.  The  mapping  from  the  physical  to  the  computational  plane  is  recomputed  be- 
cause the  Euler  solver  requires  a  different  mesh  than  that  used  by  the  potential 
solvers.  Details  of  this  new  transformation  will  be  provided  below. 

5.  The  solution  of  the  Euler  equations  is  computed,  using  the  FC  potential  solu- 
tion to  initialize. 

The  chord  of  the  airfoil  to  be  studied  is  scaled  to  one.  The  nose  is  put  at  the 
origin  and  the  trailing  edge  at  the  point  (1,0)  in  the  Cartesian  plane.  The  region 
exterior  to  the  wing  is  mapped  conformally  to  the  interior  of  the  unit  circle  of  the 
computational  plane.  The  derivative  of  the  conformal  mapping  function  is  of  the 
form 


<'-)  l=^i-r^e.pff:^.c-). 


JC  c= 


1=0 


where  z  =  x  -\-  iy,  (^  =  r  e'*,  e  is  the  (positive)  included  angle  at  the  traihng  edge, 
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and  the  Cj's  are  complex  coefficients  which  depend  on  the  geometry  of  the  airfoil. 
In  practice,  the  number  of  coefficients  computed  is  half  the  number  of  angular  mesh 
points.  The  potential  solvers  require  only  the  conformal  mapping  to  transform 
coordinates  because  the  outer  boundary  of  the  flow  domain  extends  to  infinity. 

The  Euler  solver,  however,  requires  a  mesh  which  extends  a  finite  distance 
into  the  far  field  because  it  requires  that  the  grid  boxes  of  the  physical  domain 
have  aspect  ratios  close  to  one  —  a  requirement  which  cannot  be  satisfied  by  a 
conformal  map  near  the  outer  boundary.  Consequently,  the  Euler  solver  requires  a 
second  transformation  to  rearrange  the  radicd  coordinates  of  the  grid  points  so  that 
one  may  control  the  radial  extent  and  the  aspect  ratio  of  the  far  field  grid  boxes. 
This  additional  transformation  reorders  the  radial  coordinates  of  the  grid  points  in 
the  computational  domain  so  that  the  points  nearest  the  center  of  the  unit  circle 
are  grouped  closely  together  and  thus  are  mapped  via  the  function  (3.15)  into  a 
mesh  in  the  physical  domain  which  exhibits  much  less  radial  stretching  at  the  outer 
boundary.  Figures  2  and  3  display  both  the  Euler  and  the  potential  meshes.  Note 
the  difference  in  aspect  ratio  of  the  grid  boxes  in  the  far  field.  Figure  4  displays 
the  grid  near  the  airfoil.  Note  that  the  conformed  mapping  derivative  includes  a 
Schwarz-Christoffel  term  to  allow  for  a  cusp  or  corner  at  the  trailing  edge,  and  also 
that  the  positive  segment  of  the  real  axis  in  the  (^-plane  is  opened  up  to  form  a  wake 
of  constant  thickness  in  the  2-plane.  This  is  taken  to  be  the  sum  of  the  upper  and 
lower  surface  displacements  6  at  the  trailing  edge. 
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For  an  O-mesh  with  160  angular  and  40  radial  points,  the  NC  potential  com- 
putation requires  only  a  few  seconds  on  a  Cray  Y-MP/832  supercomputer.  It  uses 
about  a  third  of  the  time  needed  by  the  FC  potential  routine  and  about  a  sixth  of 
the  time  needed  by  the  Euler.  The  FC  potential  solution  provides  a  good  initial- 
ization for  the  Euler  solver  because  theoretically  the  two  formulations  differ  only  in 
the  third  order  entropy  rise.  Obviously,  the  Euler  calculation  would  require  more 
time  if  it  did  not  make  use  of  the  potenticd  solution  as  an  initial  guess,  but  Refer- 
ence [13]  provides  additioned  techniques  for  increasing  the  speed  of  the  Euler  solver 
which  we  have  not  implemented.  In  any  case,  even  though  the  focus  of  this  work 
is  not  necessarily  on  speed,  it  is  clear  that  the  NC  potential  scheme  is  at  least  five 
times  faster  than  the  Euler,  and  it  requires  far  less  memory  space  because  only  the 
potential  (j)  must  be  stored,  rather  than  the  flow  vector  w  of  the  Euler  solver. 

Periodically  throughout  each  computation,  a  turbulent  boundary  layer  cor- 
rection is  performed  which  alters  the  shape  of  the  airfoil.  We  will  refer  to  this 
procedure  as  the  solid  displacement  method.  The  displacement  thickness  8  =  H*  6* 
is  computed  using  von  Karman's  equation 


(3.16)  !^  +  (a.+2-MM-^  =  ^ 

as  q   as       pq^ 


for  the  momentum  thickness  6*,  where  the  shape  factor  H*  and  the  skin  friction 
T  are  determined  from  semi-empiric2d  formulae  of  Nash  and  Macdonald.  A  more 
detailed  description  of  this  procedure  and  the  appropriate  references  may  be  found 
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in  Reference  [4].    The  effectiveness  of  the  solid  displacement  method  in  treating 
viscous  flow  will  be  discussed  in  Chapter  V. 

In  practice,  the  displacement  thickness  does  not  converge  to  a  steady  state,  even 
when  the  boundary  layer  correction  is  relajced.  Smzdl  oscillations  in  the  boundary 
layer  variables  do  not  have  a  great  effect  on  the  final  pressure  profile,  but  do  cause 
the  drag,  lift  and  angle  of  attack  to  noticeably  oscillate.  The  momentum  thickness 
6*  largely  determines  the  profile  drag,  the  component  of  the  total  drag  due  to  viscous 
effects,  and  thus  causes  most  of  the  uncertainty  in  the  drag  calculation.  The  wave 
drag,  which  accoimts  for  roughly  ten  percent  of  the  total  drag  in  a  typical  viscous 
calculation,  is  observed  to  converge  rather  well.  Also,  oscillations  in  the  angle  of 
attack  and  the  lift  make  it  difficult  to  compare  numerical  and  experimental  results. 
To  compensate  for  the  uncertainty  in  all  of  these  important  flow  parameters,  each 
is  calculated  as  a  weighted  average  of  the  last  few  computed  values. 
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CHAPTER  IV 
WAVE  DRAG  AND  THE  ENTROPY  INEQUALITY 

1.  Introduction 

The  addition  of  artificial  viscosity  terms  to  the  Euler  equations  must  not  only 
ensure  numerical  convergence,  but  must  also  model  the  physical  problem  correctly. 
Specifically,  the  additional  terms  must  guarantee  that  the  entropy  inequality  is  sat- 
isfied. Often  these  two  demands  may  be  met  without  having  to  add  any  terms  at  all, 
but  rather  by  implementing  a  clever  differencing  scheme.  One  example  is  presented 
by  Murman  and  Cole  [25],  in  which  upwind  differencing  in  the  supersonic  region  of 
the  flow  makes  use  of  the  truncation  error  to  represent  an  artificial  viscosity  which 
satisfies  an  entropy  inequality.  This  idea  is  used  by  the  nonconservative  potential 
solver  in  the  BGKM  code.  Its  origins  lie  in  attempts  to  solve  nonhnear  hyperbolic 
partial  differential  equations  numerically  with  difference  schemes  that  are  biased 
in  characteristic  directions.  Reference  [1]  provides  a  survey  of  the  many  numerical 
methods  which  incorporate  upwind  differencing.  The  present  Euler  solver,  however, 
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employs  central  differences  for  aJl  spatial  derivatives,  and  therefore  explicit  artificial 
viscosity  terms  must  be  added  to  aJl  four  equations  in  order  to  guarantee  conver- 
gence. A  check  on  the  physical  validity  of  the  additional  terms  is  to  show  that 
the  entropy  inequality  is  indeed  satisfied  by  the  modified  equations.  We  can  then 
derive  a  formula  for  the  wave  drag  based  on  this  inequality.  Because  this  formtda 
will  depend  on  the  specific  form  of  the  artificial  viscosity  terms,  it  will  provide  a 
further  test  of  their  Vcdidity. 

Both  the  NC  potential  and  the  Euler  solvers  in  the  BGKM  code  make  use 
of  the  entropy  inequality  to  measure  the  wave  drag.  This  calculation  involves  the 
summing  of  a  positive  definite  quantity  over  the  region  of  the  flow  in  which  the 
shock  is  smeared,  avoiding  the  use  of  computed  flow  variables  near  the  tail  which 
are  subject  to  great  uncertainty  when  a  boundary  layer  correction  is  included.  This 
is  not  true  of  the  standard  approach,  in  which  the  pressure  is  integrated  around  the 
airfoil.  Mertens,  Klevenhusen  and  Jakob  [22]  discuss  the  difficulties  of  this  standard 
method,  and  also  look  to  the  entropy  inequality  to  provide  a  more  natural  and  man- 
ageable treatment  of  the  wave  drag.  However,  their  procedure  involves  integrating 
the  jump  in  entropy  along  a  contour  arotmd  the  shock,  implicitly  assuming  that  the 
shock  will  be  sharp  and  relatively  isolated.  Reference  [4]  contains  examples  of  flows 
which  gave  rise  to  multiple  shocks  quite  near  to  each  other.  Locating  each  one  in 
order  to  carry  out  an  integration  would  be  extremely  cumbersome,  and  probably 
inaccurate  due  to  grid  spacing  considerations.    A  double  integration,  rather  than 
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a  contour  integration,  over  the  entire  shock  region  eliminates  the  need  to  locate 
individual  shocks,  enabling  us  to  treat  realistic  cases  with  multiple  shocks. 

2.  Model  Equation 

To  illustrate  the  advantages  of  computing  the  wave  drag  using  the  entropy 
inequality,  we  first  consider  Burgers'  equation.  A  similar  study  is  made  by  Bauer, 
Garabedian  and  Chang  [2]  for  the  potential  model.  Consider  the  partial  differential 
equation 

(4.1)  Wt  +  (uV2)x  =(j^ti;.)x, 

for  X  G  R  and  t  >  0,  with  u  — >  1  as  x  — >  — oo  and  u  — >  — 1  as  a;  ^  oo.   Equation 

(4.1)  represents  the  conservation  of  an  idealized  mass  flux  u^/2,  but  with  an  arti- 
ficial viscosity  term  added  on  the  right,  whose  coefficient  u  is  small,  positive,  and 
possibly  dependent  on  a;.  The  analogue  of  drag  is  the  jump  across  a  shock  of  the 
corresponding  momentum  flux  u^/3,  which  satisfies  the  steady  state  equation 

(4.2)  l—j     ={uuu;,)^-uul, 

which  is  derived  by  multiplying  (4.1)  by  u.  Integrating  this  last  equation  over  an 
interval  x  G  [a,  b],  we  obtain 


(4.3) 


u^ 


U  UUx 


b 

dx. 


Ja 


In  the  limit  as  i/  ^  0,  the  graph  of  u  approaches  a  broken  line  representing  a  shock, 
and  the  first  term  on  the  right  hand  side  of  (4.3)  goes  to  zero.  Moreover,  this  last 
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equation  verifies  the  idealized  entropy  inequzdity  u{b)  <  u{a)  across  a  shock.  The 
idealized  drag,  which  is  represented  by  the  jump  in  the  momentum  flux  u^/3,  can 
easily  be  calculated  through  the  integration  of  the  positive  definite  quantity  v  u\ 
which  is  negligible  everywhere  but  in  the  vicinity  of  the  shock. 

The  only  trouble  with  this  analysis  is  due  to  the  discretization.  If  we  were  to 
solve  equation  (4.1)  using  the  Jameson  Euler  scheme,  the  left  hand  side  of  equation 
(4.2)  would  be  differenced 


2  2        2Ax  3        2Ax  3  V       2Ax 

The  summation  of  this  term  over  a  shock  does  not  converge  exactly  to  the  the  left 
hand  side  of  (4.3)  because  at  a  discontinuity  the  second  term  on  the  right  does  not 
go  to  zero  as  Aa;  — >  0.  The  severity  of  this  error  will  be  discussed  in  Section  4. 

Obviously,  for  this  simple  problem  we  can  calculate  directly  the  idealized  drag 
from  the  values  of  u^/3  at  the  boundary  of  the  interval  [0, 1].  However,  for  the  actual 
airfoil  problem  we  cannot  accturately  calculate  the  wave  drag,  which  is  proportional 
to  the  jump  in  entropy  across  a  shock,  from  the  values  of  S  at  the  boundary  of  the 
computational  domain,  i.e.  in  the  far  field,  because  of  inaccuracies  due  to  the  mesh. 
Nor  can  we  measure  the  jump  in  entropy  by  computing  its  value  on  both  sides  of  the 
shock  because  shocks  in  general  are  smeared  over  several  mesh  points,  are  difficult 
to  locate  exactly,  are  not  necessarily  isolated,  and  are  known  to  cause  oscillations 
in  the  flow  variables  when  higher  order  schemes  are  employed.    The  method  of 
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summing  the  contributions  to  the  wave  drag  over  a  region  of  the  flow,  rather  than 
differencing  uncertain  values  of  pressure  or  entropy,  is  intended  to  circumvent  these 
difficulties. 

3.  Derivation  of  the  Wave  Drag  Formula 

As  in  the  model  problem  introduced  in  Section  2,  we  seek  to  use  the  artificial 
viscosity  terms  which  have  been  added  to  the  Euler  equations  as  a  means  of  com- 
puting wave  drag.  For  the  computation  below,  we  neglect  the  fourth  order  artificial 
viscosity  terms  because  they  are  of  third  order  in  the  grid  spacing  and  contribute 
only  a  small  amount  to  the  drag  integral  to  be  discussed  presently.  Consequences 
of  this  assumption  will  be  discussed  in  Section  4. 

Consider  the  modified  Euler  equations: 

(4.4a)  pt  +  {pu)^  +  {pv)y  =  V  •  i/Vp 

(4.4b)  {pu)t  +  {pu^  +  p)^  +  {puv)y  =  V  •  uV{pu) 

(4.4c)  {pv)t  +  (puv);,  +  (pv^  +p)y  =  V  •  uV{pv) 

(4.4d)  {pE)t  +  {puH),  +  {pvH)y  -  V  .  uV{pH), 

where  the  artificial  viscosity  operator  is  defined  as 

~  dx  \         dx  J       dy  \        dy  J 
whose  form  with  respect  to  the  transformed  Euler  equations  (3.2)  has  been  described 
in  Chapter  III.  Notice  that  the  artificial  viscosity  operator  V  •  i^V  is  apphed  to  pH 
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rather  than  pE  in  equation  (4.4d).  This  is  done  so  that  the  total  enthalpy  H  remains 
a  constant  of  integration,  even  with  the  addition  of  artificial  viscosity  terms  to  the 
equations.  The  following  derivation  is  only  slightly  altered  if  the  operator  is  applied 
to  pE  instead  of  pH,  and  we  will  present  both  forms  of  the  resulting  wave  drag 
formula. 

Equations  (4.4)  may  be  combined  to  form  a  fifth  equation  for  the  entropy 

(4.5)     {pS)t  +  (puS),  +  ipvS)y   =    ^  {  V  .  uV{pH)  -  {E-TS)V-  uVp 

-  w  (V  •  z^V(/9?/)- u  V-i/Vp)  -  v{V  ■uV{pv)-vV -uVp)  |. 

Notice  that  the  left  hand  side  of  this  equation  is  in  conservative  form.  The  right 
hand  side,  however,  may  be  separated  into  a  positive  definite  term  and  a  divergence 
term,  and  the  equation  may  be  rewritten  in  steady  state  as 


(4.6)  V  •  pSq 


CvP 


w 


WD 


+  V..{5V,+  ^  +  ^}. 


The  positive  definite  wave  drag  norm  is  defined  as 
(4.7) 


w 


WD 


|Vu||2  +  WVvW"   +  7e 


Vp 


7  +  1  /Vp    Vp'     ^ 


1      \  P       P 


Vp 


where  the  inner  product  of  two  vectors  a,  b  G  R^  is  defined  (a,  b)  =  a"^  i/  b  and  the 
corresponding  norm  is  ||a|p  =  (a,  a).  As  stated  before,  an  equation  similar  to  (4.6) 
may  be  derived  if  the  artificial  viscosity  is  appHed  to  pE  rather  than  pH  in  (4.4d). 
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Its  wave  drag  norm  is  slightly  different,  but  is  still  positive  definite: 


w 


WD 


iVulP  +  IIVvl 


+  ^'U 


Vp 


p 


7  \  P      P 


Vp 


Since  the  scheme  we  have  chosen  apphes  the  artificial  viscosity  to  pH,  we  will  have 
no  further  use  for  this  form  of  the  wave  drag  norm. 

Equation  (4.6)  may  be  integrated  over  a  region  of  the  flow  domain  D,  and 
following  an  application  of  the  divergence  theorem,  we  obtain 


(4.8)  (i  pSq-ndl  =    11^ 

an  D 


w 


WD 


dx  dy 


+    *  1/  {  SVp  +   ■■■  ]  -ndl 

BD 


To  prove  that  the  modified  Euler  equations  satisfy  the  entropy  inequality,  we  take 
D  to  be  the  arbitrarily  narrow  strip  of  width  e  surrounding  a  smeared  shock.  As 
1/  — >  0,  the  second  integral  on  the  right  hand  side  of  (4.8)  disappears  since  the 
integrand  is  bounded  on  each  side  of  the  shock,  and  we  are  left  with 


(4.9) 


f  pN  [S]  dl=   ff^  \MwD  dx  dy 

D 


where  the  integration  on  the  left  is  along  the  discontinuity,  N  is  the  normal  com- 
ponent of  the  velocity  at  the  shock,  and  the  double  integral  of  the  right  hand  side 
is  a  bounded  positive  quantity. 

Now  we  will  use  equation  (4.6)  again  to  derive  our  formula  for  computing  the 
wave  drag.  Taking  D  to  be  the  entire  flow  domain,  the  left  hand  side  of  (4.6)  reduces 
to  an  integration  at  infinity  due  to  the  no  flux  boundary  condition,  q  •  n  =  0.  The 
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second  integral  on  the  right  hand  side  may  be  neglected  since  v  is  small  and  the 
normal  derivatives  of  the  relevant  flow  variables  remain  bounded.  We  are  left  with 
the  equation 

(4.10)  (k  pSq^-n  dl  =    I      -^  \\vf\\]y^  dx  dy, 

oo  D 

which  provides  a  convenient  formula  for  calculating  the  wave  drag,  which  we  describe 
next. 

The  wave  drag,  which  is  the  component  of  the  total  drag  due  to  the  force  of 
pressure  on  the  airfoil  in  the  direction  of  the  free  stream  flow,  is  given  in  nondimen- 
sional  form  by  the  momentum  flux  integral 

(4.11)  CwD  =  ;-    (p{puq-n  +  pn^)  dl, 

Pooqio    J 

where  n^  is  the  x-component  of  the  outward-facing  normal  along  any  simple,  closed 
contour  enclosing  the  airfoil.  Let  the  free  stream  conditions  be  Uoo  =  1,  v^o  =  0, 
Poo  =  1>  and  Poo  =  7-^TO-  ^*  ^^  assumed  that  p  and  v  recover  their  free  stream 
values  in  the  narrow  strip  of  particle  paths  extending  from  the  back  of  the  shock  to 
infinity,  but  that  p  and  u  do  not.  Equation  (4.11)  may  then  be  rewritten  as 


(4.12)  CwD  =  :r    i>  pu{u  -  Uoo) 

pooqio   J 


dy. 


Since  p  is  assumed  to  return  to  its  free  stream  value  of  1  in  this  narrow  strip  at 
infinity,  u  may  be  written  as  a  function  of  the  density  p  alone 


(4.13)  u  =  Jl+^^^^         ^ 


7-1    VPoo       P 
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through  Bernoulli's  law  (2.7).    Likewise,  p  may  be  written  as  a  function  of  the 
entropy  S  alone 


(4.14)  P  =  poo  exp  (  -^ 

7C« 


through  the  equation  of  state  (2.2).  Substituting  (4.14)  into  (4.13),  we  obtain  the 
relationship  between  the  flow  velocity  u  and  the  entropy  S  at  infinity 


(4.15)  n  =  Jl  +  27eoo   |l  -  exp  (^-^^)  }• 

To  first  order  in  AS,  this  yields  the  relation 


(4.16)  Uooiu-Uoo)  =  — —{S-Soo)- 


Now  the  above  drag  formula  (4.12)  may  be  rewritten  in  terms  of  the  entropy  gain 
across  the  shock 

(4.18)  CwD  -      ^^r        ipu{S-Soo)dy. 

PooqloCv    J 
oo 

This  formula  is  presented  by  Oswatitsch  [28].  For  computational  purposes,  we 
choose  to  express  the  entropy  flux  integral  of  (4.18)  in  terms  of  the  double  integral 
of  (4.8)  computed  over  the  region  of  the  flow  domain  exhibiting  large  gradients  in 
the  flow  variables,  namely  the  neighborhood  of  the  shock.  Thus,  we  have  finally 

(4.19)  CwD  =  -^  If  -  M^wD  dx  dy. 


oo 
D 
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The  first  order  approximation  of  equation  (4.15)  leads  to,  at  worst,  about  a  3% 
error  in  the  wave  drag.  For  moderately  sized  shocks,  this  amounts  to  less  than  one 
count  of  drag  which  is  tolerable,  given  the  fact  that  the  boundary  layer  correction 
causes  far  greater  uncertainty  in  drag  prediction.  It  is  certainly  possible  to  rewrite 
equation  (4.6)  in  terms  of  a  function  F{S)  of  the  entropy  rather  than  of  S  itself,  and 
then  rederive  our  drag  formula  exactly.  However,  this  would  unduly  compUcate  the 
computation,  and  it  would  require  that  we  directly  calculate  the  jump  in  entropy 
across  the  shock  —  something  we  have  sought  to  avoid. 

4.  Implementation 

In  practice,  the  wave  drag  coefficient  is  relatively  easy  to  compute.  The  differ- 
encing of  the  flow  variables  is  dictated  by  the  way  in  which  the  artificial  viscosity 
terms  of  the  modified  Euler  equations  yield  a  summation  of  positive  definite  terms 
when  added  together  point  by  point.  We  may  thus  convert  equation  (4.19)  into 
finite  difference  form 


(4.20)    C^o  =  ^J:J:^^   {^U^tur  +  ^^IdDUr  +  •  .  .  }  Ar  A^. 


2ec^ 

,2      . 

'      j 
Note  that  DgU  and  £)+u  represent  the  forward  differences  of  u  in  each  coordinate 

direction,  the  artificial  viscosity  coefficient  u  was  defined  in  Chapter  III,  and  Ar 

and  A^  are  the  computationzd  grid  spacings. 

In  order  to  verify  that  the  formvda  (4.20)  is  scaled  properly,  the  solutions  to  the 

modified  Euler  equations  were  computed  for  the  NACA  0012  airfoil  on  meshes  of 
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various  sizes.  The  convergence  study  is  presented  in  Figure  5.  For  free  stream  Mach 
number  M  =  .78  and  angle  of  attack  a  =  0,  the  mesh  is  refined  progressively  until 
its  dimensions  are  640  by  160.  We  observe  that  the  wave  drag  coefficient  converges 
to  a  value  of  .0016,  or  16  counts  of  drag.  This  value  is  lower  than  that  of  the  FC 
potential  solver,  which  predicts  19.4  counts  using  the  standard  pressure  integration 

(4.21)  CwD  =  ^^    <p  pdy 

around  the  airfoil.  This  value  is  in  turn  lower  than  the  Euler  value  that  is  predicted 
from  the  pressure  integral  (4.21),  which  appears  to  be  heading  towards  about  22.5 
counts. 

Because  the  Euler  and  the  FC  potential  solutions  are  virtually  identical  at 
this  Mach  number,  we  expect  their  wave  drag  coefficients  to  be  the  same  as  well. 
However,  as  seen  in  the  Section  2,  the  Euler  coefficient  from  the  entropy  formula 
(4.20)  is  lower  due  to  the  discretization  error  associated  with  the  summation  of 
the  positive  definite  terms.  The  wave  drag  from  the  pressure  integral,  in  general, 
appears  to  converge  to  the  potential  value,  or  at  least  to  a  value  closer  to  it  than 
that  from  the  entropy  formula.  However,  for  the  grid  normally  used,  a  160  by  40 
0-mesh,  the  entropy  formula  usually  predicts  a  value  closer  to  the  potential.  For 
a  moderately  sized  mesh,  uncertainty  in  the  pressure  values  are  introduced  in  the 
wake  due  to  radial  stretching  and  the  higher  viscous  coefficients  required  to  enforce 
the  Kutta-Joukowski  condition.    The  wave  drag  norm  defined  by  equation  (4.7) 
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avoids  the  use  of  the  pressure  values  near  the  nose  and  tail  of  the  airfoil,  whereas 
the  pressure  integral  is  essentially  a  measure  of  their  difference. 

The  problems  with  the  pressure  integration  are  magnified  when  computations 
are  carried  out  with  a  boundary  layer  correction.  Here  we  compute  the  pressure 
on  the  boundary  of  the  displaced  airfoil,  and  thus  the  integration  around  the  air- 
foil (4.21)  must  be  replaced  by  the  corresponding  momentum  flux  integral  (4.11) 
over  a  contour  consisting  of  the  edge  of  the  boundary  layer  and  the  wake  extending 
to  infinity.  This  integration  is  unrehable  due  to  the  errors  caused  by  radial  grid 
stretching,  and  it  can  hardly  be  expected  to  provide  an  accurate  wave  drag  coeffi- 
cient. In  fact,  in  a  subsonic  viscous  computation  for  the  Korn  airfoil,  the  pressure 
integral  predicted  20  counts  of  drag,  whereas  the  entropy  formula  predicted  the 
correct  value  of  0  counts. 

The  entropy  formula  holds  even  more  promise  for  three-dimensional  compu- 
tations, because  the  wave  drag  due  to  a  3D  flow  is  an  even  smaller  component  of 
the  total  drag  due  to  the  presence  of  induced  drag,  yet  isolating  it  is  essential  for 
measuring  the  severity  of  the  shock  wave  atop  an  airfoil.  All  of  the  difficulties  in 
computing  the  wave  drag  through  a  standard  pressure  integration  which  we  have 
discussed  for  the  two-dimensional  computation  are  magnified  in  a  three-dimensional 
computation,  making  it  necessary  to  use  a  more  reliable  method  to  isolate  the  drag 
due  to  shocks.  Also,  there  is  no  reason  why  we  cannot  apply  this  new  method  to  su- 
personic or  hypersonic  flows,  where  again  the  use  of  a  standard  pressure  integration 
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is  unreliable. 

5.  Numerical  Results 

Unfortunately,  it  is  not  possible  to  test  our  new  wave  drag  formula  directly 
against  experimental  data  because  wave  drag  is  only  a  small  component  of  the  total 
drag  force  on  an  airfoil.  However,  we  may  compare  the  wave  drag  computed  by 
the  three  solvers  for  2D  inviscid  flows  past  the  NACA  0012  airfoil  at  zero  angle 
of  attack  over  a  range  of  transonic  Mach  numbers.  Figure  6  displays  the  increase 
in  the  wave  drag  with  the  free  stream  Mach  number  Moo  for  each  of  the  models 
considered.  The  NC  potential  scheme  employs  its  version  of  the  entropy  inequality 
in  calculating  the  wave  drag,  the  FC  potential  scheme  employs  the  pressure  integral 
(4.21),  and  the  Euler  scheme  employs  both  methods  for  comparison.  For  the  Euler 
solutions,  the  wave  drag  coefficients  from  the  pressure  integral  are  generally  higher 
than  the  coefficients  from  the  entropy  formula.  For  the  mesh  size  used,  it  appears 
that  the  entropy  formula  provides  a  better  estimate  of  drag  in  the  realistic  range  of 
the  Mach  spectrum. 

Figure  7  displays  the  increase  in  the  total  drag  for  viscous  computations  around 
the  Korn  airfoil  which  will  be  described  in  Chapter  V.  The  total  drag  is  the  sum  of 
the  wave  drag  and  the  drag  due  to  skin  friction.  We  compare  the  drag  coefficients 
Cd  for  the  three  numerical  solvers  with  the  experimental  drag  coefficient.  The 
Euler  solver  employs  the  entropy  formulation  to  compute  the  wave  drag,  and  this 
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leads  to  a  reasonably  good  estimate  of  the  total  drag.  Of  course,  the  component  of 
the  drag  due  to  skin  friction  dominates,  but  the  wave  drag  is  of  primary  concern  to 
the  designer  of  the  airfoil  who  seeks  to  minimize  the  drag  due  to  shocks.  For  this 
reason,  the  present  formulation  of  computing  the  wave  drag  is  a  reliable  tool  with 
which  to  compare  various  airfoil  designs. 
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CHAPTER  V 
SUPERCRITICAL  AIRFOIL  COMPUTATIONS 

1.  Introduction 

Figures  6  and  7  illustrate  that  the  wave  drag  on  a  supercritical  airfoil  rises 
quickly  with  the  free  stream  Mach  number.  In  an  actual  flow,  however,  the  boundary 
layer  separates  when  the  wave  drag  increases  above  50  counts,  and  the  model  then 
breaJts  down.  It  is  for  this  reason  that  the  isentropic  assumption  of  the  potential 
model  is  sufficient  for  applications  to  airfoils.  Since  the  entropy  rise  across  a  shock 
is  of  third  order  in  the  shock  strength,  and  since  the  shock  strength  is  bounded 
from  above  due  to  physical  considerations,  the  entropy  jump  may  be  neglected  in 
numerical  simulations  at  realistic  operating  conditions. 

In  this  chapter  we  shall  compare  the  three  models  of  flow  over  a  supercritical 
airfoil:  the  nonconservative  (NC)  potential,  the  fully  conservative  (FC)  potential, 
and  the  Euler.  Included  are  inviscid  flows  over  the  NACA  0012  airfoil  and  viscous 
flows  over  the  Kom  airfoil  for  various  free  stream  Mach  numbers  and  angles  of 
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attack.  To  model  viscous  effects,  we  implement  a  boundary  layer  correction  using 
the  solid  displacement  method  which  has  been  described  in  Chapter  III.  Schmidt, 
Jameson,  and  Whitfeld  [31]  employ  another  method,  known  as  the  surface  source 
method,  which  we  will  describe  below.  To  our  knowledge,  no  attempt  has  yet  been 
made  to  couple  the  Jameson  Euler  scheme  with  the  solid  displacement  boundary 
layer  correction. 

Due  to  inaccuracies  introduced  by  the  mesh  used  by  the  Euler  solver  (see 
Chapter  III,  Section  4),  all  of  the  Euler  plots  to  follow  exhibit  a  discrepancy  in 
the  coefficient  Cp  which  shifts  the  pressure  profile  to  a  slightly  higher  value.  This 
shift  comes  as  the  result  of  radial  stretching  of  grid  cells  and  does  not  affect  the 
main  thrust  of  the  study,  which  is  primarily  focused  on  the  discrepancy  in  chordwise 
shock  location  of  the  three  methods. 

2.  The  NACA  0012  Airfoil  and  Shock  Location 

Figure  8  presents  the  results  of  inviscid  simulations  carried  out  over  the  sym- 
metrical NACA  0012  airfoil  at  zero  angle  of  attack  for  Mach  numbers  ranging  from 
a  subsonic  value  M  —  .74  to  the  highly  supercritical  case  M  =  .85  .  Notice  that 
the  subsonic  solutions  predicted  by  the  three  models  are  essentially  identical.  They 
also  remain  close  during  the  onset  of  transonic  flow  at  the  Mach  numbers  .76  and 
.78  .  Supercritical  airfoils  under  normal  operating  conditions,  that  is,  airfoils  whose 
boundary  layers  are  unseparated,  typically  exhibit  shocks  of  comparable  strength  to 
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those  displayed  in  these  three  figures.  The  purpose  of  these  inviscid  calculations  is 
to  compare  how  each  model  predicts  shock  location  and  wave  drag.  In  Chapter  III 
we  presented  the  wave  drag  results,  and  we  now  consider  the  issue  of  shock  location. 
Let  us  postpone  for  the  moment  consideration  of  the  NC  potential  solution 
and  focus  our  attention  instead  on  the  conservative  methods,  which  differ  only  in 
regard  to  the  conservation  of  momentum  and  entropy  across  a  shock.  Since  the 
entropy  jump  across  a  shock  is  of  third  order  in  the  shock  strength,  it  is  reasonable 
to  expect  that  the  difference  in  shock  location  between  the  FC  potential  solution 
and  the  Euler  solution  is  small  in  realistic  cases.  This  is  indeed  the  case,  judging 
from  the  succession  of  pressure  plots  in  Figure  8.  Figure  9  presents  the  relationship 
between  the  jump  in  entropy  and  the  difference  in  shock  location;  both  measured 
at  the  airfoil  surface.  We  represent  the  entropy  gain  by  the  quantity 

where  A(S)  is  the  coefficient  used  to  relate  pressure  and  density  in  the  equation  of 
state  (2.2),  and  the  subscripts  '2'  and  '1'  denote  values  of  S  on  the  airfoil  boundary 
behind  and  in  front  of  the  shock,  respectively.  Even  though  these  measurements  are 
subject  to  considerable  error  from  the  oscillations  and  smearing  which  occur  at  the 
shock,  we  are  fairly  confident  that  the  shift  in  shock  location  is  of  the  same  order 
in  the  shock  strength  as  the  gain  in  entropy. 
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3.  The  Korn  Airfoil  with  Boundary  Layer  Correction 

Figures  10  through  15  display  the  pressure  profiles  for  flows  past  the  Korn 
airfoil,  computed  with  a  boundary  layer  correction.  No  attempt  has  been  made 
to  match  the  angles  of  attack  with  those  of  the  experimental  tests  with  which  the 
calculations  are  compared  [15].  Instead,  the  free  stream  Mach  number  Moo  and  the 
coefficient  of  lift  Cl  were  chosen  as  the  proper  parameters  with  which  to  match 
the  numerical  to  the  wind  tunnel  data.  We  do  this  because  the  experimental  angle 
of  attack  is  subject  to  significant  error.  The  supercritical  simulations  carried  out 
earlier  by  Bauer,  Garabedian,  Korn,  and  Jameson  [4]  match  lift  coefficient  rather 
than  angle  of  attack,  and  the  success  of  their  results  validates  our  procedure. 

The  inviscid  simulations  of  Section  2  show,  for  the  conservative  methods,  a 
steep  recovery  in  the  pressure  behind  the  shock  which  is  known  as  the  Oswatitsch- 
Zierep  singtilarity.  It  has  been  observed  that  severe  gradients  of  this  kind  cannot 
be  sustained  in  the  presence  of  a  boundary  layer  [9],  and  consequently,  the  pressure 
jump  in  a  viscous  flow  tends  to  be  smaller  than  that  predicted  by  the  Rankine- 
Hugoniot  shock  conditions.  Shock  wave-boundary  layer  interaction  is  the  key  to 
understanding  this  phenomenon. 

The  von  Karman  equation,  used  in  the  BGKM  code  to  compute  the  bound- 
ary layer,  predicts  that  a  shock  will  cause  the  displacement  thickness  6  to  increase 
rapidly.  This  eff"ect  is  handled  differently  by  the  two  conventional  methods  of  in- 
corporating a  turbulent  boundary  layer  correction  into  flow  computations.    The 
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first,  known  as  the  solid  displacement  method,  adds  a  thickness  S  to  the  airfoil  and 
computes  the  inviscid  flow  past  the  altered  body.  The  second  method,  based  on 
a  surface  source  idea  of  Lighthill  [18],  simulates  a  boundary  layer  by  changing  the 
no-flux  boundary  condition  at  the  airfoil,  allowing  mass  to  enter  the  flow  domain 
according  to  the  rule 


(5.2)  pq-n=  —  (p5(5), 


where  n  is  the  outward  normal  to  the  surface  and  5  is  a  measure  of  arclength. 

The  BGKM  code  employs  the  solid  displacement  method,  and  in  cases  where 
the  shock  falls  farther  back  on  the  airfoil,  the  NC  potential  solution  is  in  better 
agreement  with  the  experimental  data  than  is  either  of  the  conservative  methods. 
The  FC  potential  and  the  Euler  solutions  predict  shocks  which  are  larger  and  farther 
back  than  the  experimental  shock  and  the  shock  predicted  by  the  nonconservative 
model.  Perhaps  the  reason  for  this  is  that  the  slight  increase  in  the  displacement 
thickness  at  the  shock,  which  is  significant  in  a  real  flow,  is  essentially  not  detected 
in  a  numerical  integration  of  the  von  Karman  equation  which  tends  to  smooth  out 
small  disturbances.  It  appears  that  the  mass  production  associated  with  the  NC 
potential  equation  at  a  shock  [9]  thus  models  shock  wave-boundary  layer  interaction 
quite  well.  The  neglect  of  the  jump  in  entropy  by  the  potential  methods  is  of 
secondary  importance  in  these  cases. 

The  NC  potential  model  yields  a  less  accurate  solution,  however,  when  the 
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shock  develops  farther  forward  on  the  airfoil,  as  in  Figure  11(a).  An  examination  of 
the  computed  values  of  the  displacement  thickness  in  Figure  11(b)  shows  that  the 
boundary  layer  is  relatively  thin  and  does  not  undergo  a  dramatic  thickening  near 
the  weak  shock.  Thus  the  shock  wave-boundary  layer  interaction  is  less  important 
when  shocks  are  formed  where  the  boundary  layer  is  relatively  thin.  Consequently, 
the  conservative  methods  here  provide  a  more  appropriate  model  of  shock  structure. 
The  Euler  option,  in  particualar,  satisfies  the  Rankine-Hugoniot  relations  described 
in  Chapter  II,  but  yields  a  solution  which  is  more  or  less  equivalent  to  the  FC 
potential  solution  —  the  isentropic  assimiption  is  again  of  minor  consequence. 
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CHAPTER  VI 


VORTEX  DYNAMICS 


1.  Introduction 

The  preceding  chapter  verifies  what  has  been  known  for  some  time  —  that 
the  potential  model  provides  the  aeronautical  engineer  with  an  efficient  tool  with 
which  to  design  and  evaluate  supercriticcd  wing  sections.  However,  the  fact  that 
the  Euler  equations  admit  rotational  solutions  makes  them  of  interest  for  other 
appHcations  as  well.  In  this  chapter  we  consider  rotational,  inviscid  flow  around 
airfoils  and  various  blunt  bodies.  Again,  by  inviscid  we  mean  that  boundary  layer 
effects  are  not  considered  —  the  modified  Euler  equations  obviously  are  equipped 
with  viscous  terms.  Of  great  concern  to  our  study  is  vortex  tracking,  which  is  usually 
accomplished  numerically  by  computing  the  paths  traced  out  by  point  vortices  with 
no  dissipative  mechanisms  present.  The  Euler  solver  we  have  been  considering  is 
able  to  track  vortices  reasonably  well  before  its  diff'usive  terms  weaken  and  wash 
them  downstream. 
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The  focus  here  is  not  to  model  the  most  difficult  physical  phenomena,  such  as 
turbulence,  but  rather  to  provide  insight  into  the  artificial  viscosity  mechanism  of 
the  scheme.  For  the  discussion  to  follow,  we  will  make  reference  to  the  modified 
Euler  equations  written  in  idealized  form 

(6.1)  wt  +  f(w)^  +  g(w)j,  =  V  ■  1/2 Vw  -  V  •  1/4V  Aw. 

The  vectors  are  all  defined  in  Chapter  III.  We  have  seen  that  the  second  order  term  of 
the  artificial  viscosity  is  needed  to  satisfy  the  entropy  inequality  in  the  neighborhood 
of  a  shock,  where  the  coefficient  1^2  is  of  first  order  in  the  grid  spacing.  Now  we 
shall  see  that  the  fourth  order  term,  whose  coefficient  V4  is  of  third  order,  is  not 
only  needed  to  stabilize  the  computation  in  the  far  field,  but  also  to  control  vortices 
and  to  impose  the  Kutta-Joukowski  condition  at  the  traihng  edge. 

The  flows  presented  in  this  chapter  are  also  of  interest  because  they  may  be 
computed  without  modifying  the  airfoil  code  or  refining  the  mesh.  By  simply  chang- 
ing the  magnitude  of  the  artificial  viscosity  control  parameters  K2  and  K4,  we  may 
compute  the  flow  around  a  wide  variety  of  objects  ranging  from  the  circular  cylinder 
to  a  flat  plate  perpendicular  to  the  free  stream.  We  will  even  demonstrate  the  flex- 
ibility of  the  scheme  by  simulating  such  well  known  phenomena  as  the  von  Karman 
vortex  street  and  the  Helmholtz  instability. 

We  should  mention  again  that  these  flow  solutions  are  advanced  in  time  using 
the  local  CFL  condition.  Because  the  vorticity  is  usually  found  near  the  body,  the 
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vaxiation  in  the  time  step  does  not  affect  the  solution  too  much.  However,  as  the 
vortices  extend  farther  from  the  solid  boundary,  the  non-uniformity  of  the  mesh 
plays  a  greater  role  in  their  movement, 

2.  The  Circular  Cylinder  and  the  Problem  of  Uniqueness 

The  Euler  equations  admit  an  infinite  number  of  rotational  solutions  in  which 
there  occur  closed  streamlines  which  do  not  reach  the  flow  boundary  or  extend  to 
infinity.  Even  though  the  Euler  equations  which  we  solve  numerically  have  had 
artificial  viscosity  terms  added,  the  question  of  uniqueness  must  still  be  addressed. 
Since  the  viscosity  coefficients  1/2  and  U4  are  both  positive,  the  artificial  viscosity 
terms  of  (6.1)  are  clearly  dissipative.  Vorticity  tends  to  be  diff"used  and  washed 
downstream,  causing  converged  solutions  to  be  irrotational,  or  nearly  so.  Even 
when  shocks  or  sharp  edges  serve  as  sources  of  vorticity,  the  artificial  viscosity 
terms  often  limit  their  effect  on  the  overall  solution.  In  Chapter  V  we  have  shown 
that  the  vorticity  due  to  weak  shocks  on  airfoils  is  effectively  diffused  by  the  artificial 
viscosity,  rendering  the  solutions  more  or  less  irrotational.  Now  we  will  study  this 
phenomena  for  blunt  body  flows. 

With  no  mechanism  to  create  vorticity,  i.e.  a  boundary  layer,  shock,  or  sharp 
edge,  a  flow  which  is  initially  isentropic  and  irrotational  remains  that  way.  The 
modified  Euler  equations  in  this  case  can  certainly  be  expected  to  yield  the  po- 
tential flow  solution,  and  indeed  they  do.  If,  on  the  other  hand,  the  computation 
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is  initialized  with  a  rotational  flow,  we  only  see  the  emergence  of  the  irrotational 
solution  after  a  sufficient  number  of  iterations  have  been  performed  to  diffuse  the 
vorticity  and  wash  it  downstream. 

We  have  computed  solutions  to  the  modified  Euler  equations  for  the  flow  past 
a  circular  cyhnder,  initializing  the  computation  with  the  incompressible  flow  which 
includes  two  symmetric  point  vortices  in  the  wake.  In  terms  of  the  complex  coor- 
dinate 2  =  X  +  iy,  the  complex  potential  of  this  flow  is  given  by 


a 


2 


(6.2)  <l)  +  ixl;  =  z+ i  riog 


{z  -  zo){z  -  a^/zo) 


z  (z  -  zo){z  -  a^/zo)' 

where  a  is  the  radius  of  the  cyUnder  and  zq  and  zq  are  the  locations  of  point 
vortices  of  strength  ±r.  We  have  selected  these  parameters  to  conform  with  a 
stable  stationary  solution  [23],  but  we  have  smoothed  out  the  singularity  in  the 
velocity  at  the  points  zo  and  zq.  Because  we  are  at  the  moment  interested  in 
subsonic,  rotational  flows,  we  have  imposed  a  free  stream  Mach  number  M  =  .2  so 
that  we  may  turn  off  the  second  order  difference  terms  of  the  artificial  viscosity  by 
setting  the  parameter  ^2  =  0.  Of  course,  in  the  absence  of  a  shock  the  second  order 
coefficient  would  have  been  small,  but  we  want  to  study  here  how  the  fourth  order 
terms  alone  handle  the  vorticity. 

Figure  16  displays  two  flow  patterns  computed  with  different  values  of  the 
parameter  K4.  A  relatively  low  Vcdue  is  used  in  16(a)  and  a  higher  value  in  16(b). 
Notice  in  16(a)  that  a  Helmholtz  instability  has  arisen  because  the  backward  flow 
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at  the  boundary  of  the  cyUnder  slips  under  the  forward  flow,  but  in  16(b)  the 
flow  separates  from  the  boundary  at  a  weU  defined  point.  Obviously,  the  larger 
coefficient  of  the  fourth  order  terms  prevents  the  backwash  from  occurring.  In 
either  event,  both  computations  will  eventually  converge  to  the  potential  solution 
after  the  vorticity  is  diffused  and  washed  downstream. 

We  may  also  initialize  the  flow  with  the  vortices  in  front  of  the  cyhnder.  Figure 
17  displays  the  evolution  in  time  of  the  vortex  pair.  At  first  their  respective  strengths 
are  large  enough  to  propel  them  upstream,  but  eventually  they  are  weakened  and 
washed  back  downstream.  Since  the  entropy  S  is  convected  by  the  flow,  we  see  the 
emergence  of  vortices  behind  the  cyUnder  as  those  in  front  disappear.  This  example 
is  of  interest  to  designers  of  turbine  or  helicopter  blades,  who  may  not  take  the  free 
stream  flow  to  be  isentropic  since  each  blade  follows  in  the  wake  of  another. 

If  the  free  stream  Mach  number  is  large  enough  to  cause  a  shock  to  develop,  then 
as  predicted  by  the  Rankine-Hugoniot  shock  conditions  and  by  Crocco's  equation 


(6.2)  ^-qxu;  =  TVS-VH, 

at 


an  initially  irrotational  flow  becomes  rotational.  For  the  computational  model  we 
are  considering,  the  shock  is  the  principal  mechanism  which  introduces  vorticity 
into  an  initially  irrotational  flow,  for  there  is  no  sharp  edge  or  boundary  layer  to 
create  vorticity.  The  vortices  are  diffused  by  the  numerical  viscosity  and  convected 
downstream  by  the  flow.  This  should  result  in  a  continuous  transport  of  vorticity 
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from  the  body  by  way  of  the  wake.  Presented  in  Figure  18,  however,  is  a  series  of 
stills  of  a  cyUnder  flow  in  which  the  mechanism  for  vorticity  production,  the  shock, 
and  the  mechanism  for  vorticity  depletion,  the  artiflcicd  viscosity,  are  out  of  step. 
Thus  the  vorticity  created  by  the  shock  causes  two  large  vortices  to  form  behind 
the  cyhnder  which  so  dramatically  aJter  the  apparent  shape  of  the  body  that  the 
shock  loses  its  strength  and  disappears.  Then  with  no  source  of  vorticity,  the  large 
vortices  are  diffused  by  the  artificial  viscosity  and  swept  downstream,  and  the  shock 
reappears  as  the  body  regains  its  cyhndrical  shape.  This  sequence  of  events  repeats 
itself  over  and  over  so  that  the  flow  never  attains  a  steady  state.  The  example 
calls  into  question  the  abihty  of  the  modified  Euler  equations  to  simulate  correctly 
the  behavior  of  a  flow  which  becomes  rotational  due  to  a  shock.  The  simulations 
of  Chapter  V  succeed  in  modeUng  inviscid  flows  over  airfoils  because  K4  is  large 
enough  to  dissipate  the  vorticity  rapidly,  before  the  body  shape  is  altered  by  large 
vortices. 

An  interesting  consequence  of  this  last  example  is  the  von  Karman  vortex 
street.  If  we  continue  to  iterate  the  flow  pictured  in  Figure  18,  eventually  round-off 
error  causes  the  solution  to  lose  symmetry.  Figure  19  presents  a  few  stills  in  which 
the  vortices  are  now  shed  at  intermittent  time  intervals.  Unfortunately,  the  viscous 
dissipation  washes  out  the  vortices  so  that  it  is  not  possible  to  view  the  entire  train 
of  alternating  swirls. 
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3.  The  Kutta-Joukowski  Condition 

Kelvin's  theorem  asserts  that  an  isentropic,  inviscid  flow  originating  in  a  free 
stream  will  remain  irrotational.  Thus,  in  a  reaJ  time-dependent  flow,  the  circulation 
established  around  an  airfoil  due  to  the  Kutta-Joukowski  condition  is  counterbal- 
anced by  a  small  vortex  of  equal  magnitude  rotating  in  the  opposite  direction,  which 
is  shed  off'  the  trailing  edge  as  the  airfoil  begins  its  movement  from  rest.  The  po- 
tential model  predicts  the  resultant  steady  state  flow  by  explicitly  enforcing  the 
Kutta-Joukowski  condition  through  the  boundary  condition  (2.9)  and  through  the 
requirement  that  the  velocity  be  finite  at  the  traihng  edge.  These  conditions  estab- 
lish a  circulation  at  infinity  to  match  the  circulation  established  around  the  airfoil. 
Because  the  circulation  is  prescribed  at  infinity  through  a  far  field  boundary  con- 
dition, the  potential  solvers  are  extremely  fast  and  robust,  allowing  for  either  the 
coefficient  of  lift  or  the  angle  of  attack  to  be  fixed  at  the  outset  of  the  computation. 

The  time-dependent  Euler  solver,  however,  includes  no  such  explicit  mechanism 
for  enforcing  the  Kutta-Joukowski  condition.  It  has  often  been  asserted  that  the 
computed  Euler  solution  automatically  satisfies  this  requirement  [31].  We  have 
found  that  this  is  not  really  true,  for  the  Euler  scheme  relies  on  the  fourth  order 
difference  terms  of  the  artificial  viscosity  to  suppress  motion  of  the  fluid  around 
the  trailing  edge  of  the  airfoil,  and  thereby  produce  a  solution  satisfying  the  Kutta- 
Joukowski  condition.  In  order  to  demonstrate  this,  an  inviscid  computation  was 
conducted  for  the  Korn  airfoil  at  zero  angle  of  attack,  starting  with  a  uniform  flow 
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as  the  initial  condition.  Figure  20(a)  displays  the  flow  field  which  is  established 
at  the  traihng  edge  after  only  19  iterations.  This  flow  pattern  resembles  those  of 
actual  physical  experiments  in  which  a  small  vortex  forms  just  above  the  trailing 
edge  of  the  airfoil  and  is  then  shed  off  as  circulation  is  established  around  the 
airfoil.  Unfortunately,  the  computation  fails  during  the  next  iteration  because  of 
its  inability  to  treat  the  cavitation  created  by  a  rapid  increase  in  the  flow  velocity 
around  the  tail  of  the  airfoil. 

There  are  two  remedies  for  this  problem.  The  first  is  to  increase  the  viscosity 
parameter  K4,  which  controls  the  fourth  order  terms,  so  that  it  becomes  large  enough 
to  suppress  the  irregularity  created  at  the  traihng  edge.  For  the  above  numerical 
experiment,  k,4  =  .005.  If  we  increase  it  to  .015,  we  see  in  Figure  20(b)  that 
the  vortex  is  suppressed,  and  in  Figure  20(c)  finally  disappears  when  the  value  is 
increased  to  .025.  The  other  remedy  to  the  instability  of  Figure  20(a)  is  to  initialize 
the  Euler  solution  with  the  potential  solution.  This  eliminates  comphcations  arising 
from  the  original  vortex  shedding,  but  does  not  eliminate  the  difficulty  altogether 
because  changes  in  the  lift  continually  alter  the  circulation  around  the  airfoil.  This 
problem  is  remedied  again  by  an  increase  in  the  fourth  order  artificial  viscosity 
parameter.  In  either  case,  one  must  choose  a  proper  value  of  K4  to  impose  the 
Kutta-Joukowski  condition. 

We  see  from  the  above  discussion  how  the  Kutta  condition  is  enforced  at  the 
wing  by  the  artificial  viscosity,  but  there  remains  the  issue  of  the  far  field.    The 
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Euler  solver  makes  no  use  of  the  circulation  around  the  airfoil  in  prescribing  the  far 
field  boundary  conditions.  For  this  reason  it  is  extremely  difficult  to  prescribe  the 
lift,  rather  than  the  angle  of  attack,  in  a  flow  computation.  The  time  lag  required 
for  information  at  the  airfoil  surface  to  propagate  to  the  outer  boundary  makes  it 
hard  to  iterate  the  angle  of  attack  for  a  fixed  lift.  Thomas  and  Salas  [32]  discuss 
the  problems  associated  with  the  circulation  and  the  far  field  boundary  conditions, 
concluding  that  if  the  outer  boundary  of  the  mesh  for  an  Euler  computation  is  not 
sufficiently  far  from  the  airfoil,  say  at  least  10  chord  lengths,  it  becomes  necessary  to 
impose  a  circulation  at  the  outer  boundary  based  on  the  small  disturbance  equation. 
This  is  more  or  less  what  is  accomplished  by  the  potential  solvers.  However,  the 
mesh  used  in  the  computations  described  in  Chapter  V  extends  roughly  25  chord 
lengths  from  the  airfoil,  far  enough  to  eliminate  the  need  to  impose  a  circulation  at 
the  outer  boundary. 

4.  The  Stalled  Wing 

The  Euler  flow  around  a  stalled  wing  sheds  further  light  on  the  discussion  of 
Sections  2  and  3.  Figure  21  presents  a  typical  flow  pattern  for  the  NACA  0012 
airfoil  at  an  angle  of  attack  of  15  degrees.  As  with  extreme  cases  of  flow  around  a 
circular  cyfinder,  we  have  not  been  able  to  compute  a  unique  steady  state  solution 
for  the  airfoil  in  buffet  because  of  the  interaction  of  the  mechanisms  for  vorticity 
production  and  depletion.  Here  vorticity  is  created  by  the  shock  at  the  nose  of  the 
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airfoil  and  also  by  the  sharp  trailing  edge. 

Let  us  again  consider  the  Kutta-Joukowski  condition.  Figure  22(a)  shows  the 
velocity  vectors  at  the  traihng  edge  for  a  flow  simulation  with  a  low  value  of  K4. 
Notice  that  the  fluid  does  not  leave  the  body  in  the  direction  tangent  to  the  lower 
surface  due  to  the  interference  of  a  vortex  which  has  formed  just  off^  the  upper 
surface.  The  next  plot,  Figure  22(b),  shows  the  effect  on  the  flow  of  doubling  the 
parameter  K4.  Now  the  fluid  trails  off  the  body  in  the  right  direction,  and  the 
vortex  formed  at  the  interface  of  the  converging  flows  does  not  interfere  with  the 
flow  emanating  from  the  lower  surface.  The  fourth  order  artificial  viscosity  terms 
essentially  prevent  vortices  from  forming  on  small  scales  by  smoothing  out  large 
gradients  in  the  flow  variables.  Vortices  form  on  large  scales,  however,  as  can  be 
seen  in  Figure  21.  These  large  scale  vortices  do  not  cause  any  problems  for  the  flow 
solver,  even  at  the  tail,  as  long  as  the  artificial  viscosity  parameters  are  suflRciently 
large. 

5.  Persistent  Rotational  Flows 

We  conclude  with  two  cases  where  the  solution  does  not  converge  to  the  po- 
tential solution,  but  rather  remains  rotational  due  to  the  presence  of  a  source  of 
vorticity. 

As  a  preliminary  study  of  the  applicability  of  the  scheme  to  model  three- 
dimensional  rotational  flows,  we  present  the  flow  formed  by  inserting  a  vertical 
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plate  in  a  uniform  2D  flow.  This  gives  us  insight  into  how  similar  codes  handle  the 
vortex  sheet  shed  from  the  wing  tips  of  an  airplane  in  a  3D  calculation.  Figure  23 
displays  the  flow  field  for  a  typical  case.  The  flat  plate  is  in  fact  an  ellipse  with  5% 
aspect  ratio  that  is  tipped  up  90  degrees.  We  use  a  thin  ellipse  rather  than  a  flat 
plate  to  avoid  making  any  modification  to  the  existing  algorithm.  The  flow  pattern 
displayed  in  Figure  23  is  precisely  what  one  would  hope  to  see.  Because  the  flow 
goes  around  a  relatively  sharp  edge,  the  artificial  viscosity  cannot  wash  away  the 
large  vortices  shed  by  the  plate,  and  we  arrive  at  the  steady  state  flow  because  there 
is  a  constant  source  of  vorticity. 

As  we  have  seen  in  Section  2,  the  Euler  solver  is  capable  of  treating  the  inter- 
face between  oppositely  directed  streams,  simulating  the  Helmholtz  instability.  An 
example  of  this  phenomenon  is  presented  in  Figure  24.  We  initialize  the  flow  past 
a  thin  ellipse  with  free  stream  values  v  =  0,  p  =  1,  p  =  l/M^,  Moo  =  -5,  and 

(1       y  <yo 
-1       y  >  yo, 

where  yo  is  selected  so  that  the  interface  is  far  enough  from  the  airfoil  to  avoid  early 
interference  with  the  solid  boundary.  Figure  24  displays  the  entire  mesh  which 
extends  about  25  chord  lengths  from  the  airfoil,  thus  the  ellipse  is  shrouded  by 
the  fine  mesh  surrounding  it.  The  specific  body  shape  is  not  of  interest  here,  but 
rather  the  abihty  of  the  scheme  to  compute  the  flow  on  a  non-uniform  mesh.  At  a 
suitable  intermediate  time,  a  virtually  periodic  train  of  vortices  is  observed  of  the 
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kind  seen  in  other  calculations  in  which  periodic  initial  data  have  been  imposed 
over  a  periodic  mesh.  The  flow  pattern  clearly  shows  how  the  artificial  viscosity 
smooths  out  large  gradients  in  the  velocity.  Eventually,  the  string  of  vortices  will 
be  diffused,  and  the  resulting  flow  will  consist  of  large  sweeping  streams  from  the 
incoming  to  the  outgoing  parts  of  the  outer  boundary,  with  vorticity  being  created 
at  the  solid  surface  boundary  at  the  origin  of  the  flow  field. 

We  have  included  these  last  two  examples  in  our  study  merely  to  stress  the 
usefulness  of  the  Euler  solver  in  handling  rotational  flow.  Even  though  the  potential 
model  is  preferable  in  treating  the  airfoil  problem,  it  cannot  be  used  if  rotatioued 
effects  are  important.  For  example,  in  a  3D  aircraft  computation  the  induced  drag 
due  to  vortex  shedding  at  the  wing  tips  is  of  major  concern.  Whether  or  not  the 
Euler  solver  can  simulate  this  phenomenon  remains  to  be  seen.  The  studies  in  this 
chapter  suggest  that  the  artificial  viscosity  mechanism  of  the  Jameson  scheme  would 
certainly  be  a  valuable  tool  for  these  applications. 
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Figure  1.  Typical  analysis  run  of  the  NC  potential  solver. 
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Figure  2.  Grid  used  by  the  Euler  solver. 
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Figure  3.  Grid  used  by  the  potential  solvers. 
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Figure  4.  Grid  around  the  displaced  Korn  airfoil. 
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Figure  5.  Convergence  study  for  wave  drag  calculations. 
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Figure  6.  Drag  rise  for  the  NACA  0012  airfoil. 
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Figure  7.  Drag  rise  for  the  Korn  airfoil. 
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Figure  8.  (a)  Subsonic  inviscid  flow  computation  for  the  NACA  0012  airfoil. 
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Figure  8.  (b)  Inviscid  flow  computation  for  the  NACA  0012  airfoil. 
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Figure  8.  (c)  Inviscid  flow  computation  for  the  NACA  0012  airfoil. 


73 


-1  .2 


.6   - 


1  .2 


1  .B 


CP 


NACA12  AIRFOIL 

—  NC  POTENTIAL 

—  FC  POTENTL^L 

—  EULER 


M*N=  160*40  NO  VISCOSITY 

M  =  .800     ALP=  0.00  CL=0.000     CD  =  .0063 

M=.800     ALP=  0.00  CL^O.OOO     CD  =  .0080 

M=.800     ALP=  0.00  CL=0.000     CD  =  .0062 


Figure  8.  (d)  Inviscid  flow  computation  for  the  NACA  0012  airfoil. 
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Figure  8.  (c)  Inviscid  flow  computation  for  the  NACA  0012  airfoil. 
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Figure  8.  (f)  Inviscid  flow  computation  for  the  NACA  0012  airfoil. 
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Figure  9.  Shift  in  shock  location  for  Euler  vs.  FC  potential. 
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Figure  10.  Viscous  flow  computation  —  subsonic  case. 
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Figure  11.  (a)  Viscous  flow  computation  —  forward  shock  case. 
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Figure  11.  (b)  Upper  surface  displacement  thickness  6. 
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Figure  12.  Viscous  flow  computation  —  slightly  off  design. 
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Figure  13.  Viscous  flow  computation  —  high  Mach  number  case. 
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Figure  14.  Viscous  flow  computation  —  high  Mach  number  case. 
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Figure  15.  Viscous  flow  computation  —  Korn  airfoil  in  buffet. 
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Figure  16.  Circvdar  cylinder:  (a)  low  K4.  (b)  high  K4. 
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Figure  17.  (a)  Initial  solution,  (b)  Vortex  pair  moves  upstream. 
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Figure  17.  (c)  Vortex  pair  is  diffused,  (d)  Vortices  reappear  behind. 
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Figure  18.  Shock-vortex  interaction. 
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Figure  19.  Onset  of  the  von  Karman  vortex  street. 
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Figure  20.  Kutta  condition:  (a)  K4  =  .005  (b)  K4  =  .015  (c)  K4  =  -025. 
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Figure  21.  NACA  0012  airfoil  tipped  up  15  degrees. 
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Figure  22.  (a)  Divergence  at  the  trailing  edge,  (b)  Stability. 
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Figure  23.  Flow  past  a  vertical  plate. 
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Figure  24.  The  Helmholtz  instability. 
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